Many X Variables

Guidance Counselors

Next, let’s recreate some visualizations of these regression lines.

plot(df$hs_gpa, df$sat_sum, las = 1, pch = 19,
     col = scales::alpha(ifelse(df$sex == 1, "tomato", "dodgerblue"), 0.4),
     ylab = "SAT Percentile Sum",
     xlab = "High School GPA")
abline(r2, col = "tomato", lwd = 3)
abline(r3, col = "dodgerblue", lwd = 3)
abline(r1, col = "black", lwd = 2, lty = 2)
legend("bottomright", bty = "n", lty = c(2, 1, 1),
       col = c("black", "tomato", "dodgerblue"),
       legend = c("All", "Female", "Male"), lwd = 2)

Splitting our sample like this is convenient because it allows for a lot of flexibility – each group gets their own parameters. However, remember: the more parameters you need to estimate, the more difficult inference becomes.¹

If we wanted to reduce the number of parameters in our model (which is currently four: two intercepts and two slopes), we could simply not differentiate by gender. By doing this, we are sacrificing a bit of bias for reductions in variance. Then, we’d only have the two parameters (intercept and slope) to estimate. Is there a compromise we can make, and maybe estimate just three parameters?

Guidance Counselors

The short answer is: yes! Let’s explore this a bit.

In the two-model situation, we have two intercepts and two slopes. In order to reduce the number of parameters, we would need to make an assumption about the setting we’re working with. We could either have male and female students share a slope or an intercept. If they share a slope, we would only have to estimate the single slope parameter and two intercept parameters. After looking at the two regression lines, they seem to be pretty parallel, which means that the same-slope assumption seems reasonable.

Guidance Counselors

Before we estimate the model, let’s run through some of the algebra of this model. Suppose we are working with a male student, and they have a 3.0 GPA. What do we expect their SAT to be?

What about a female student with the same GPA?

The only difference in expected SAT scores is \(\widehat{\beta}_1\), and this is independent of \(\text{GPA}\). In other words, we would interpret \(\beta_1\) as an intercept shifter/modifier, or the average difference between female and male scores.

Let’s estimate the model and plot the resulting lines.

Guidance Counselors

# create binary indicator
df$female <- ifelse(df$sex == 1, 1, 0)
r4 <- lm(sat_sum ~ female + hs_gpa, df)
summary(r4)

Call:
lm(formula = sat_sum ~ female + hs_gpa, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.642  -8.401  -0.083   8.540  34.198 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  59.9776     2.4687  24.295   <2e-16 ***
female        6.8978     0.7927   8.702   <2e-16 ***
hs_gpa       12.4561     0.7335  16.982   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.39 on 996 degrees of freedom
Multiple R-squared:  0.248, Adjusted R-squared:  0.2465 
F-statistic: 164.2 on 2 and 996 DF,  p-value: < 2.2e-16

1. The intercept is the expected sum of SAT percentiles when hs_gpa and female are equal to zero. Therefore, the interpretation for this coefficient is now the expected sum of SAT percentiles for a male student with a GPA of zero.
2. Keeping hs_gpa at zero, the expected sum of SAT percentiles for female students is 59.98 + 6.9 = 66.88.
3. Increasing hs_gpa has the same effect for both groups given our initial assumption. A one unit increase in hs_gpa leads to an increase in SAT of 12.46

Guidance Counselors

To allow for differing slopes by group, we can specify the following functional form:

So, in other words, the intercepts would be exactly the same, which matches our assumption.

Guidance Counselors

What if these two students had 3.0 GPAs? How would their expected SAT scores change?

This makes the difference between the two groups equal to \(3\widehat{\beta}_2\), or \(\text{GPA} \times \widehat{\beta}_2\). If \(\widehat{\beta}_2\) is positive (negative), the difference between the groups would widen (shrink) as \(\text{GPA}\) increased.

Guidance Counselors

To estimate this regression in R, we will write the equation as follows: sat_sum ~ hs_gpa + hs_gpa:female. The colon (:) tells R to multiply the two variables in the model like we have written above.

r5 <- lm(sat_sum ~ hs_gpa + hs_gpa:female, df)
summary(r5)

plot(df$hs_gpa, df$sat_sum, las = 1, pch = 19,
     col = scales::alpha(ifelse(df$sex == 1, "tomato", "dodgerblue"), 0.4),
     ylab = "SAT Percentile Sum",
     xlab = "High School GPA")
abline(a = coef(r5)[1], b = coef(r5)[2], col = "dodgerblue", lwd = 3)
abline(a = coef(r5)[1], b = sum(coef(r5)[2:3]), col = "tomato", lwd = 3)
legend("bottomright", bty = "n", lty = c(1, 1),
       col = c("tomato", "dodgerblue"),
       legend = c("Female", "Male"), lwd = 2)

Call:
lm(formula = sat_sum ~ hs_gpa + hs_gpa:female, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-43.884  -8.284  -0.255   8.600  34.832 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    63.9518     2.3804  26.866  < 2e-16 ***
hs_gpa         11.3040     0.7288  15.511  < 2e-16 ***
hs_gpa:female   2.0289     0.2446   8.294 3.54e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.43 on 996 degrees of freedom
Multiple R-squared:  0.2431,    Adjusted R-squared:  0.2415 
F-statistic: 159.9 on 2 and 996 DF,  p-value: < 2.2e-16

Here, the two lines certainly diverge as hs_gpa increases, which matches up with the regression coefficient (2.03). To interpret this coefficient: female student SAT scores grow two points more than male student SAT scores when increasing GPA by one.

The coefficient for the product (multiplication) of two variables in a regression is called an interaction term.

Guidance Counselors

Econometricians like to use dummy variables and interaction terms because they allow us to consider the setting’s inherent heterogeneity in our models. However, another convenient thing about them is that they allow us to estimate differences.

For example, what if we wanted to know: do female students preform better on SAT exams than male students? Before this module, we might collect SAT scores for female and male students then do a difference-in-means hypothesis test (from Module 2). If we found a statistically significant difference, however, we could not be sure that the difference was in fact due to gender but rather due to the female students in the sample just being smarter overall. By using regression, we can partial out the effect of intelligence (via GPA¹) and then re-calculate the difference in means. In fact, our female dummy variable from r4 is exactly that: a difference in means test. The difference we found was 6.9, and it was statistically significant even after controlling for hs_gpa.

Guidance Counselors

Now that we’ve estimated these two models with three parameters each, we can combine them into one giant model. This model will look like the following:

Coefficient Interpretation:

1. \(\beta_0\) is the expected SAT for male students with a GPA of 0.
2. \(\beta_1\) is the difference in the expected SAT for female students relative to male students with a GPA of 0.
3. \(\beta_2\) is the change in expected SAT for male students if they increase their GPA by 1.0
4. \(\beta_3\) is the difference in the change in expected SAT for female students relative to male students when they increase their GPA by 1.0.

Housing Price Models

Of course we love pretending to be guidance counselors, but we really love being real estate agents. For (maybe?) the last time, let’s read in our ames housing data (data; documentation). We are going to keep all of the variables we did last time plus one more: BsmtFin.Type.1.

BsmtFin Type 1  (Ordinal): Rating of basement finished area

       GLQ  Good Living Quarters
       ALQ  Average Living Quarters
       BLQ  Below Average Living Quarters   
       Rec  Average Rec Room
       LwQ  Low Quality
       Unf  Unfinshed
       NA   No Basement

ames <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/openintro/ames.csv")
keep_columns <- c("price", "area", "Bedroom.AbvGr", "Year.Built", "Overall.Cond", "BsmtFin.Type.1")
ames <- ames[,keep_columns]
colnames(ames) <- c("price", "sqft", "bedrooms", "yr_built", "condition", "bsmnt")
ames$age <- 2011 - ames$yr_built

Housing Price Models

Let’s take a peak at the bsmnt variable to see what all it contains. To do this, we are going to use table(), but add in useNA = "always" just in case there are any missing data elements in the column.

table(ames$bsmnt, useNA = "always")

      ALQ  BLQ  GLQ  LwQ  Rec  Unf <NA> 
   1  429  269  859  154  288  851   79 

According to this tabulation, there appears to be 79 observations with no input, and 1 observation with a blank input. According to the definition above, NA values mean “no basement”.¹ Let’s change the data from NA to "No Basement".

ames$bsmnt <- ifelse(is.na(ames$bsmnt) | ames$bsmnt == "", "No Basement", ames$bsmnt)

Housing Price Models

If we were real estate agents, it would be really nice to have some way to quantify this basement variable so we could use it in our model. Otherwise, we might end up like Zillow with a potentially huge omitted variable issue. To start, let’s just create a dummy variable equal to 1 if the property has a basement and a zero otherwise.

ames$has_bsmnt <- ifelse(ames$bsmnt == "No Basement", 0, 1)

Similar to how we controlled for gender differences, which is a categorical variable, we can include this variable into our housing price model. For simplicity, we are going to just use square footage as our other explanatory variable. Students are encouraged to explore the inclusion of other variables by themselves.

Housing Price Models

Estimating our model:

reg1 <- lm(log(price) ~ log(sqft) + has_bsmnt, ames)
reg1

Call:
lm(formula = log(price) ~ log(sqft) + has_bsmnt, data = ames)

Coefficients:
(Intercept)    log(sqft)    has_bsmnt  
     5.1596       0.8931       0.3881  

This model suggests that a 1% increase in square footage increases expected sale price by 0.89%. In addition, having a basement increases price by 38.8%. This is helpful, but there are many different “ratings” of basement quality given in bsmnt that we are ignoring. The lowest quality is “Unf”, or unfinished. Let’s make a dummy for this type of basement and include it in the model.

Housing Price Models

This interpretation is a bit clunky because it requires us to add up two coefficients to get our total estimate. Rather, to make interpretation easier, we can redefine has_bsmnt to exclude unfinished basements:

ames$has_bsmnt <- ifelse(ames$bsmnt %in% c("No Basement", "Unf"), 0, 1)
reg2 <- lm(log(price) ~ log(sqft) + has_bsmnt + unfinished, ames)
reg2

Call:
lm(formula = log(price) ~ log(sqft) + has_bsmnt + unfinished, 
    data = ames)

Coefficients:
(Intercept)    log(sqft)    has_bsmnt   unfinished  
     5.0663       0.9062       0.4363       0.2683  

This interpretation is much easier. An unfinished basement, relative to no basement, is worth a 26.8% increase. A non-unfinished basement, relative to no basement, is worth a 43.6% increase. Hopefully you can start to see a pattern. If we wanted to estimate the effect of some other basement quality level, we would have to again redefine has_bsmnt and create a new dummy variable. If we take this logic to the finish line, we’d end up \(N-1\) dummy variables where \(N\) is the number of possible categories. We subtract one because one group must be the reference group.

Housing Price Models

I’m sure you’ve heard this one before: work smarter, not harder. In programming, they say: A good programmer is a lazy programmer. While one of these is earnest and the other sarcastic, they mean similar things. In our case, making a bunch of dummy variables can be a lot of work! Luckily, R has a lazy efficient solution for us.

To create many dummies all at once, we can use the factor() command. Moreover, we can even select which factor becomes our reference group! To do this in R, we would do the following:

# I am putting a "_" to differentiate between the original variable.
ames$bsmnt_ <- as.factor(ames$bsmnt)
head(ames$bsmnt_)
cat("\n")
ames$bsmnt_ <- relevel(ames$bsmnt_, ref = "No Basement")
head(ames$bsmnt_)

[1] BLQ Rec ALQ ALQ GLQ GLQ
Levels: ALQ BLQ GLQ LwQ No Basement Rec Unf

[1] BLQ Rec ALQ ALQ GLQ GLQ
Levels: No Basement ALQ BLQ GLQ LwQ Rec Unf

Housing Price Models

The model spits out a bunch of numbers for bsmnt_, but it’s important to remember that they are all interpreted as having that “type” of basement relative to not having a basement, controlling for the property’s age, size and number of bedrooms. Moreover, we can plot the coefficients and show that they’re generally increasing in the order that we’d expect.

BsmtFin Type 1  (Ordinal): Rating of basement finished area

       GLQ  Good Living Quarters
       ALQ  Average Living Quarters
       BLQ  Below Average Living Quarters   
       Rec  Average Rec Room
       LwQ  Low Quality
       Unf  Unfinshed
       NA   No Basement

It is important to point out the assumptions of this model. We are assuming that changing the square footage (or age or bedrooms) of a home impacts sale price independently of the type of basement. Of course, we could estimate a fully interacted model, but the output would be huge and uninterpretable.

Fixed Effects

One of the most common ways to use dummy variables is for fixed effects. A common use case for fixed effects is when you are working with a time-by-unit panel. For example: a state-by-year panel.

By including state dummies, or fixed effects, we are removing / controlling for any and all constant (time-invariant) factors. As an example, if the outcome variable in our state-by-year panel is housing prices, we can expect California to have high prices and West Virginia to have low prices. Why? For a variety of reasons: weather and other local amenities, economic factors, etc. To make the states more comparable to one another, we can remove the effect of the idiosyncratic differences intrinsic to each state.

We can also include year fixed effects, we can remove any and all factors that impact each state at the same time. For example, the 2007 financial crisis reduced housing prices across the country. When we include dummies for each year (except one, the reference group), we can remove the influence of these events on our outcome.

Fixed Effects

To demonstrate fixed effects, we are going to compare per-capita personal income and housing price indices for different counties/cities in Hampton Roads over time. Let’s take a look at the data files we have.

hr <- read.csv("https://alexcardazzi.github.io/econ311/data/hr_fips.csv")
hr # hr is "hampton roads"

                                    name  fips
1                      Gloucester County 51073
2                   Isle of Wight County 51093
3  James City County + Williamsburg City 51931
4                         Mathews County 51115
5                            York County 51199
6                        Chesapeake City 51550
7                           Hampton City 51650
8                      Newport News City 51700
9                           Norfolk City 51710
10                       Portsmouth City 51740
11                          Suffolk City 51800
12                   Virginia Beach City 51810

Fixed Effects

To begin merging, we should check out which data set has the most “coverage”. For example, the pi dataset contains the range of years 1969, 2021 whereas hp contains only 1975, 2022. This can also be seen graphically:

par(mfrow = c(1, 2))
plot(table(pi$year), las = 1,
     main = "Personal Income",
     xlab = "Year", ylab = "Frequency")
plot(table(hp$year), las = 1,
     main = "Housing Price Index",
     xlab = "Year", ylab = "Frequency")
par(mfrow = c(1, 1))

These figures suggest we should merge into pi, or else we would potentially be dropping data. It’s always a good idea to have too much than too little.

Fixed Effects

Now let’s use the same code to pull the FIPS names from hr.

m <- match(pi$fips, hr$fips)
pi$name <- hr$name[m]
View(pi[,c("value", "fips", "year", "name")])

Fixed Effects

There appears to be a pretty strong, positive relationship between these variables. Now let’s take a look at what some regressions tell us about the strength of the relationships. First, we are going to estimate a model where housing price index is the outcome and personal income is the explanatory variable. Next, we are going to include FIPS and year fixed effects.

coef(lm(hp ~ value000, pi))[2]; cat("\n")
coef(lm(hp ~ value000 + as.factor(fips) + as.factor(year), pi))[2]

value000 
3.698684 

 value000 
0.4945878 

Why are these two coefficients so different?

Linear Probability Model

So far in this module, we have learned about binary explanatory variables. However, what if your outcome variable is binary? We can still use OLS for this!

Let’s consider some data on police treatment of individuals arrested in Tortonto for simple possession of small quantities of marijuana (data, documentation). In these data, we can observe whether the arrestee was released, their race, age, sex, employment status, citizen status and previous police interactions.

Let’s build a model to explain whether the arrestee was released given some of these controls.

df <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/carData/Arrests.csv")
df$released <- ifelse(df$released == "Yes", 1, 0)
df$black <- ifelse(df$colour == "Black", 1, 0)
df$female <- ifelse(df$sex == "Female", 1, 0)
df$citizen <- ifelse(df$citizen == "Yes", 1, 0)
df$employed <- ifelse(df$employed == "Yes", 1, 0)
summary(r_lpm <- lm(released ~ female + citizen + employed + age + black + checks, df))

Call:
lm(formula = released ~ female + citizen + employed + age + black + 
    checks, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.98272  0.03631  0.09415  0.19252  0.60548 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.7405159  0.0238808  31.009  < 2e-16 ***
female      -0.0034265  0.0179731  -0.191    0.849    
citizen      0.0892898  0.0143901   6.205 5.89e-10 ***
employed     0.1236378  0.0125926   9.818  < 2e-16 ***
age          0.0004879  0.0006049   0.807    0.420    
black       -0.0573499  0.0119889  -4.784 1.77e-06 ***
checks      -0.0499783  0.0033967 -14.714  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3581 on 5219 degrees of freedom
Multiple R-squared:  0.09551,   Adjusted R-squared:  0.09447 
F-statistic: 91.85 on 6 and 5219 DF,  p-value: < 2.2e-16

Linear Probability Model

For variables with few unique values, like checks, we can transform them into factors even though they are not technically categorical. This can help us capture non-linear effects. Below are estimated coefficients from a linear probability model with check fixed effects. In addition is a plot of the coefficients with their 95% confidence intervals to visualize the results.

r_ck <- lm(released ~ black + female + citizen + employed + age + as.factor(checks), df)
summary(r_ck)

Call:
lm(formula = released ~ black + female + citizen + employed + 
    age + as.factor(checks), data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.96115  0.05103  0.09370  0.18417  0.57696 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         0.7313119  0.0240721  30.380  < 2e-16 ***
black              -0.0605972  0.0119997  -5.050 4.57e-07 ***
female              0.0001721  0.0179659   0.010   0.9924    
citizen             0.0884361  0.0144127   6.136 9.09e-10 ***
employed            0.1221697  0.0125839   9.708  < 2e-16 ***
age                 0.0003205  0.0006079   0.527   0.5981    
as.factor(checks)1 -0.0155505  0.0148940  -1.044   0.2965    
as.factor(checks)2 -0.0452317  0.0154331  -2.931   0.0034 ** 
as.factor(checks)3 -0.1315396  0.0148573  -8.854  < 2e-16 ***
as.factor(checks)4 -0.2152762  0.0168947 -12.742  < 2e-16 ***
as.factor(checks)5 -0.2534383  0.0330511  -7.668 2.07e-14 ***
as.factor(checks)6 -0.1208529  0.1200205  -1.007   0.3140    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3575 on 5214 degrees of freedom
Multiple R-squared:  0.09935,   Adjusted R-squared:  0.09745 
F-statistic: 52.29 on 11 and 5214 DF,  p-value: < 2.2e-16

Reporting Coefficients

When reporting regression models, people are usually uninterested in fixed effects. Rather, people like to include table rows that indicate which fixed effects are included in which models. As an example, let’s consider the housing price and personal income regressions.

fe1 <- lm(hp ~ value000, pi)
fe2 <- lm(hp ~ value000 + as.factor(fips) + as.factor(year), pi)

regz <- list(`Housing Price Index` = fe1,
             `Housing Price Index` = fe2)
coefz <- c("value000" = "Personal Income (in Thousands)")
gofz <- c("nobs", "r.squared")

rowz <- as.data.frame(matrix(c('Fixed Effects', "No", "Yes",
                               "Something Else", "A", "B"),
                             byrow = T, nrow = 2)) # 2 rows
attr(rowz, 'position') <- c(3, 4) # rows 3 and 4

modelsummary(regz,
             title = "Effect of Personal Income on Housing Price Index",
             estimate = "{estimate}{stars}",
             coef_map = coefz,
             gof_map = gofz,
             add_rows = rowz)